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1 Introduction 


The Higgs-like boson discovered at LHC [1, 2] initiated great efforts to pin it down as a particle 
responsible for electroweak symmetry breaking. Although the measured properties of this new bo¬ 
son are consistent with the expectations for the Standard Model (SM) Higgs boson [3], a variety 
of beyond-SM interpretations is still possible. One such possibility is the interpretation as a light 
state within a richer spectrum of scalar particles in the theoretically well motivated Minimal Su¬ 
persymmetric Standard Model (MSSM). Its Higgs sector consists of two complex scalar doublets 
leading to five physical Higgs bosons and three (pseudo-) Goldstone bosons. The physical mass 
eigenstates at the tree level are the neutral CP-even h, H, the CP-odd A, and the charged H ± 
bosons. Their masses can be parameterized by the charged (or CP-odd) Higgs-boson mass and 
the ratio of the two vacuum expectation values, tan/3 = v^/vx- CP-violation in the Higgs sector 
is induced by complex parameters in other sectors of the MSSM via loop corrections, leading to 
mixing between h, H , and A in the mass eigenstates [4-7]. 

Since the masses of the neutral Higgs bosons are strongly affected by loop contributions, a lot 
of work has been invested into higher-order calculations of the mass spectrum from the SUSY 
parameters, in the case of the real MSSM [8-23] as well as for the MSSM with complex parame¬ 
ters (cMSSM) [5-7, 24-29]. At one-loop order, the dominant contributions arise from the Yukawa 
sector with the large top-Yukawa coupling h t , or a t = h% / {An). The class of leading two-loop 
Yukawa-type corrections of O(a^) has recently been calculated for the case of complex parame¬ 
ters [28, 29]. Together with the full one-loop result [26] and the leading 0{atCt a ) terms [27] it is 
now available in the public program FeynHiggs [10, 20, 26, 30-32]. 

The 0(aj ) contributions in the cMSSM have been evaluated in the Feynman-diagrammatic ap¬ 
proach, extending the on-shell renormalization scheme of Ref. [26] to the two-loop level. This ensures 
that the obtained analytical results for the renormalized two-loop self-energies can consistently be 
incorporated in FeynHiggs. For implementation in FeynHiggs the whole calculation of the 0(a1) 
contributions was reworked and broken up into several working steps. Together they provide a 
semi-automated framework for a two-loop calculation in a model with non-trivial renormalization. 
Although the code is specific to the evaluation of the mentioned corrections, it can be utilized as a 
template for similar calculations. The aim of this paper is to explain in detail how the code works 
and what is done at each step, and where process-specific adjustments are necessary. 

The paper is organized as follows: Sect. 2 provides the theoretical foundation of the computa¬ 
tion of the Higgs-mass corrections in the cMSSM. The approximations used for the evaluation of 
the O(a^) contributions to the Higgs masses are explained in Sect. 3. A detailed description of the 
implementation of the 0(aq:) corrections in FeynHiggs is given in Sect. 4. 


2 Masses of the Higgs bosons in the Complex MSSM 


Since the calculation of the 0(al) corrections has already been presented in great detail else¬ 
where [28, 29] we shall recap only the parts necessary for the implementation in FeynHiggs here. 
In the basis of the tree-level Higgs states the loop-corrected mass matrices are given by [26] 
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where the E denote the renormalized self-energies. In the following we shall neglect the Goldstone 
matrix elements (greyed out above) since their numerical effect is tiny. Mixing with Goldstone 
bosons is taken into account inside the loop diagrams and for a consistent renormalization, of 
course. The physical masses are obtained as the real parts of the poles of the propagator matrix 

A(p 2 ) = -i[p 2 l-M 2 ] _1 . (2.3) 

This is tantamount to finding the zeroes of the determinant of p 2 l — M 2 , which is core functionality 
of FeynHiggs. 

For the cMSSM, the following contributions to Eq. (2.1) are available in FeynHiggs: the one-loop 
self-energies with full p 2 dependence [26], the leading momentum-independent two-loop 0(a s at ) 
self-energies [27], and now also the 0(a 2 ) [29] self-energies. The subleading two-loop 0(a s ctb), 
O(atab) parts are implemented only in the rMSSM yet, and are interpolated in the phases when 
complex parameters are chosen [26] . 

For the latter corrections the following components have to be computed: 

• The unrenormalized genuine two-loop self-energies Ej^, E^, E^, E^, E^, E^, yS^ + h _ 
at p 2 = 0 in 0{a 2 ) approximation. 

• The one-loop diagrams with insertions of one-loop counterterms. For the 0(af) corrections 
the couplings and masses of the colored sector and chargino-neutralino sector are affected by 
this subrenormalization. 

• The two-loop counterterms for these self-energies. The counterterms, the imposed renormal¬ 
ization conditions, and all renormalization constants are given in Ref. [29]. 

• The two-loop tadpoles T f [ 2 \ T^\ in 0(a 2 ) approximation, which enter the two-loop 
counterterms. 

3 Approximations in the C>(a^) contributions 

The dominant terms of the 0(a 2 ) contributions are enhanced by an additional factor m 2 . The 
following approximations are applied to the renormalized two-loop self-energies in Eq. (2.1) to yield 
only this leading part. 

3.1 Gaugeless limit 

The gauge couplings g\ and <72 are set to zero, and also the strong coupling g s is discarded. Conse¬ 
quently, also the gauge-boson masses tf^ and Mz vanish, while the weak mixing angle 9 W retains 
its original value. These choices cannot naively be substituted in the original MSSM model file 
of FeynArts [33] where the couplings are not yet expressed in terms of the top-Yukawa coupling 
ht = errit/{s/2 sp s w Mw)- Furthermore, the non-zero ratios 5M\y/M\v and SMz/Mz appear in 
the renormalization of h. t , while the gauge-boson-mass renormalization constants SMw and SMz 
themselves vanish. The technical implementation of these issues is described in Sect. 4. 

3.2 Vanishing external momentum 

The external momentum is set to zero, i.e. the integrals that appear in the two-loop self-energies and 
the renormalization constants are calculated at p 2 = 0. As an immediate consequence all appearing 
two-loop integrals are independent of p 2 , leading to vacuum diagrams that are known analytically. 
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3.3 Bottom mass equal to zero 

For the 0(a 2 ) contributions to form a supersymmetric and gauge-invariant class of two-loop cor¬ 
rections to the Higgs-boson masses, the bottom-quark mass has to be set to zero. 


3.4 Consequences of the approximations 

The approximations mentioned above modify the MSSM parameters that enter the calculation as 
follows: 


• The Higgs- and Goldstone-boson masses at the tree level become m\ = tuq = m 2 G ± = 0 and 
m 2 H =m\ = m 2 H± . 

• The Higgs-sector mixing angles a and /3 are related through a = /3 — tt/2. 

• The stop masses are derived as the eigenvalues of 


/ to? l + ro? rritX't \ 

y m t X t m? R + ml) ’ 


(3.1) 


with X f = A t — n*/tp, where the soft-SUSY-breaking parameters m 2 ^ r and A t enter. 

• Only one sbottom remains in the Feynman diagrams with m~ = ■ The sbottom mixing 

matrix is set to the unity matrix. 

• Of the neutralinos and charginos only the original Higgsinos Xg 4 , xt occur hr the Feynman 
diagrams and their degenerate masses are equal to \/i]. The mixing matrices N for the neu¬ 
tralinos and U, V for the charginos simplify to 
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• Many couplings are simplified. Technically this is taken care of by a modified version of 
the FeynArts model file for the complex MSSM, generated as described in Sect. 4. Explicit 
results can also be found in Ref. [29]. 


4 Implementation in FeynHiggs 

The original calculation of the renormalized two-loop self-energies at 0(a 2 ) [29] was overhauled 
entirely for inclusion in FeynHiggs. The new code is included in the FeynHiggs tarball and, while 
still specific to the corrections mentioned, is very adaptable and may serve as a template for similar 
calculations, which is the reason we describe it here in some detail. 

The calculation has been divided into seven working steps implemented by the scripts described 
in the following sections. Additionally, several packages that are required for specific purposes are 
accessed during each step. A makefile coordinates the entire calculation. The total running time 
is about 15-20 minutes and the final output is a highly optimized and very compact Fortran code 
with a file size of about 350 kBytes. The result can directly be moved into the FeynHiggs source 
tree. 

The advantages of this setup are twofold: First, the compartmentalization takes care of ‘con¬ 
current’ software packages (which cannot be loaded into the same Mathematica session, typically 
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because of symbol conflicts). Likewise, if an extension of this code required a different package for 
a particular step (say, tensor reduction), only this one step would need to be modified, and the 
additional package could be an arbitrary one, not just a Mathematica package, since it is invoked 
from a shell script. Second, the makefile sees to it that only the necessary parts of the calculation 
are redone in case of a change, which was a particular time-saver during the debugging phase. 
Some details of the software setup: 

• The code described below is included in FeynHiggs 2.11.0 or later, in the gen/tlsp tree. 

• The shell scripts in the ‘scripts’ subdirectory perform the actual calculation. They are 
executed directly from the command-line (i.e. not loaded within Mathematica), usually by 
the makefile. Internally they run the Mathematica Kernel as described in Ref. [34]. 

• The Mathematica programs in the ‘packages’ subdirectory provide functions for specific pur¬ 
poses and are read by the scripts. 

• As in the FeynArts model file for the MSSM [33] we choose the particle labels 

h = hO, H = HH, A = AO, H + = Hp, H~ = Hm. 

• Except for the last step (Code Generation) each script operates on a single self-energy or 
tadpole only, which is referred to on the command-line and stored under one of the following 
labels: 

self-energies: hOhO, hOHH, hOAO, HHHH, HHA0.A0A0, HmHp, 
tadpoles: hO, HH, AO. 

• Symbolic expressions are generally output to the ‘m’ directory, Fortran code to the ‘ f ’ directory. 

• Input and output files are defined in statements ‘in=. ..’ and ‘out=. . .’ in the first few lines 
of each script. For debugging, a script’s session dialog is written to the output file with 
.log.gz appended. For example, m/hOhO/6-comb. log.gz is the log of the run that made 
m/hOhO/6-comb. 

4.1 Model-file preparation, Gaugeless Limit 

Our calculation uses the FeynArts model file for the MSSM with counterterms [35]. Since the 
calculation is carried out in the gaugeless limit (cf. Sect. 3), it speeds up calculations to prepare a 
variant of the model file where this limit is taken already at the level of the Feynman rules. At this 
stage we also introduce h t (where hf = 47 t at) and substitute A t by X t + n*/tp. 

The original MSSMCT.mod model file contains both the counterterm vertices and the definitions 
of the (one-loop) renormalization constants. We moved the latter to model/MSSMCT.rcl to make 
modification of the Feynman rules easier. 

For each step we summarize usage as follows: 

• Command line: model/O-glmod 

• Main packages used: packages/Gaugeless .m, packages/XtSimplify ,m 

• Input: model/MSSMCT.mod.in 

• Output: model/MSSMCT.mod, model/MSSMCTgl.mod 
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Figure 1. The 1PI two-point two-loop topologies with propagator numbers. 


4.2 Step 1: Diagram generation 

Diagrams are generated with FeynArts [36] using the MSSMCTgl. mod model file prepared above. Cus¬ 
tomized wrapper functions for the FeynArts functions are provided by packages/FASettings ,m. 

Diagram selection in particular is simplified by filters defined in FASettings .m. Defining a 
diagram filter which works for an arbitrary topology is not entirely straightforward at the two-loop 
level but since only a few topologies contribute we have chosen a per-topology approach where 
particle choices can be made for each propagator. 

For example, the two-loop self-energy diagrams for the neutral Higgs bosons are picked with the 
following statement: 

sel [0] [S[_] -> S[_]] = { 
t [3] kk htb [6] , 
t [3] kk tb [6] , 
t [3] kk tb [6] , 
t [3] kk t [4] kk htb [5] , 
t [3] kk htb [5 I 6] , 
t [3] kk htb [5] , 
t[3] kk t[5] , 
t [5] kk ht [3 I 4] , 
t [3|4|5] && ht [3 I 4 I 5] > 

The nine elements of the list correspond to the nine two-loop topologies of Fig. 1. An entry like 
t [3] mandates that propagator //3 must carry a top or stop. Likewise, htb [516] enforces a 
Higgs/Higgsino, top/stop, or bottom/sbottom on propagators #5 and/or $=6. 
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• Command line: scripts/l-amps se CTO 

se is one of the self-energy/tadpole names given above (hOhO, HHHH, etc.) and cto is the 
counterterm order, i.e. 0 for the virtual two-loop diagrams and 1 for the one-loop plus counter¬ 
term diagrams. 

• Main packages used: packages/FASettings.m, packages/Gaugeless .m 

• Input: (none) 

• Output: m/sE/l-amps. cto (amplitude), m/sE/l-amps. cto . ps. gz (diagrams) 

4.3 Step 2: Preparation for Tensor Reduction 

In preparation for Step 3 (Tensor Reduction) several simplifications are carried out: Firstly, the 
limit p 2 —> 0 is applied. This limit is part of the approximation employed here and moreover admits 
closed-form expressions (logs and dilogs) for the two-loop integrals. Secondly, simplifications are 


done in particular on the ubiquitous sfermion mixing matrix elements 1 Uij = USf Li,jl ( i,j = 1,2). 
Products are substituted as 

UCSf Li , j] = USf Li,jl USf Li , j] *, (4. 1) 

UCSf [i, 3] = USf [z , 1] USf [i , 2] * , UCSf [3,3] = USf [1 , 1] USf [2,2] *, (4.2) 

UCSf [3,j] = USf [l,j] USf [2,j]*, UCSf [3,4] = USf [1 ,2] USf [2, 1] * (4.3) 

which makes it a lot simpler to exploit unitarity: the relation UU' = 1, for example, directly 


leads to UCSf [1,1] = UCSf [2,2] and UCSf [1,2] = UCSf [2,1]. A subtle technical aspect is that 
not all unitarity relations could even be formulated as definitions for the USf (being too ‘deep’ for 
Mathematica’s UpValues) and applying them as rules would not work reliably since Mathematica 
has no particular incentive to arrange the USf in just the way to make the match. 

After applying unitarity it is observed that almost all UCSf in the amplitude appear in the form 


of only four linear combinations: 

U2sl[x] = (UCSf [1,3] x* + UCSf [1,3]* x)/2, (4.4) 

U2s2[x] = (UCSf [1,3] x* -UCSf [1,3]* x)/(2i), (4.5) 

U2c 1 [x] = (UCSf [3 , 3] x* + UCSf [3,4]* x)/2 , (4.6) 

U2c2 [x] = (UCSf [3 , 3] x* - UCSf [3,4]* x)/(2i). (4.7) 


Substituting these combinations indeed shortens the amplitude by about 1/3. Again, this choice is 
empirical and may well be specific to the corrections computed here. All relations for substitution 
and simplification are contained in packages/U2Simplify ,m. Note in particular that, despite the 
suggestive notation, Eqs. (4.6) and (4.7) are complex quantities since the two summands are not 
conjugates of each other. 

• Command line: scripts/2-prep se cto 
se and cto as in Step 1. 

• Main packages used: packages/USfSimplify,m, packages/U2Simplify,m, 
packages/SimplificationDefinitions,m 

• Input: m/ se/ 1-amps . cto 

• Output: m/sE/2-prep. cto 

1 The USf actually carry two more indices for sfermion type and generation: USf [/ ,j, t , , which we suppress here. 
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4.4 Step 3: Tensor Reduction 

The tensor reduction of the loop integrals is performed with the packages TwoCalc [37, 38] (virtual 
two-loop diagrams, cto = 0) and FormCalc [39] (one-loop plus counterterm diagrams, cto = 1). 

One feature of this script is that, because TwoCalc and FormCalc cannot be loaded into one 
Mathematica session due to symbol conflicts, the shell script invokes not one but two different 
Mathematica sessions depending on the value of cto. 

• Command line: scripts/3-calc se cto 
se and cto as in Step 1. 

• Main packages used: TwoCalc, FormCalc, packages/FCSettings ,m, 
packages/UseSimplePackage.m 

• Input: m/sE/2-prep. cto 

• Output: m/sE/3-calc. cto 

4.5 Step 4: Simplification after Tensor Reduction 

The tensor reduction is traditionally the step that increases the number of terms in the amplitude 
most. Step 4 reduces the size of each individual amplitude before they are combined in Step 6. 

The key to decent performance comes from a ‘tag’ function 2 DiagMark [mj which the FeynArts 
wrapper functions of Step 1 (packages/FASettings .m) inserted for each diagram, where rrq are the 
masses that run in the loop. This DiagMark function induces a partitioning, i.e. the entire amplitude 
can be written as the sum of pieces, each of which is multiplied by a different DiagMark. Typically 
only few simplifications can be made across these pieces, hence one can restrict application of a 
simplification function simp to each piece, as in 

Collect [amp, _DiagMark, simp] 

This is of course much faster than applying simp to the entire expression amp. 

• Command line: scripts/4-simp se cto 
se and cto as in Step 1. 

• Main packages used: packages/Simplif icationDef initions .m, packages/FCSettings ,m 

• Input: m/sE/3-calc. cto 

• Output: m/sE/4-simp. cto 

4.6 Step 5: Calculation of the Renormalization Constants 

The one-loop Renormalization Constants (RCs) are defined in MSSMCT.rcl (see Sect. 4.1), where 
the 8My (V = W, Z) require special attention, however. As already pointed out in Sect. 3, they 
are proportional to My and would naively vanish in the gaugeless limit, while 8 My/My must 
not vanish. Technically, we accomplish this by explicitly writing out the proportionality to My\ 
dMVsql —> dMVsqlMV2*MV2, such that a possible MV2 in the denominator of the coefficient of dMVsql 
can cancel. 

For the two-loop RCs of the cMSSM no complete set of definitions is available as a model file 
yet. 3 We added the necessary ones in the approximation relevant here in model/MSSMCT.rc2. For 
better extendability we tried to preserve the mathematical notation as much as possible, e.g. the 
two-loop wavefunction renormalization is defined via (the mass counterterm is defined separately) 

2 The DiagMark function was originally developed for Refs. [40, 41], 

3 A complete list for the counterterms of the Higgs sector in the complex MSSM can be found in Ref. [29]. 
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dMHmat[Msq_, dZl_, dZ2_, dMl_] := 

1/2 (Msq.dZ2 + ConjugateTranspose[dZ2].Msq) + 

1/2 (ConjugateTranspose[dZl].dMl + ConjugateTranspose[dMl].dZl) + 

1/4 ConjugateTranspose[dZl] .Msq.dZl 

which preserves the matrix notation of Eq. (2.30) of [29]. 

The RCs are further expanded in the dimension parameter e = (D — 4) /2. For example, the 
one-loop top-mass RC, dMf 1 [3,3], is decomposed as 

dMf1[3,3] = RC[-1, dMf1[-1,3,3]] + RC[0, dMfl[0,3,3]] 

where the two terms on the r.h.s. represent the coefficients of £ _1 and £°, respectively. The loop 
integrals are expanded with the ExpandDel function and powers of e are collected with DelSeries, 
both defined in packages/ExpandDel ,m. The notation allows to easily extract the terms actually 
needed in the final result in Step 6. 

The O (e 1 ) parts of the one-loop integrals often cancel in the final result. This cancellation should 
(and in the present case has been) checked, yet it needs some effort to see the coefficients of the e 1 - 
terms actually vanish. Once this is established, however, it is useful to set the e^terms to zero manu¬ 
ally, which is faster and gives a more compact result. In packages/ExpandDel. m this is done through 
rules defined in AssumeVanishingDel. To re-enable checking, set AssumeVanishingDel = {}. 

The actual calculation of the RCs is carried out with FormCalc, for each self-energy and tadpole. 
So as not to re-compute RCs many times, they are cached in the m/rc subdirectory. Take care that 
this cache is not controlled by the makefile. 

• Command line: scripts/5-rc se 
se as in Step 1. 

• Main packages used: packages/RenConst .m, packages/ExpandDel ,m 

• Input: m/sE/4-simp.1 

• Output: m/sE/5-rc, m/rc/* 

4.7 Step 6: Combining Results 

The final analytical step for each self-energy/tadpole joins the amplitudes with the renormalization 
constants and, unless a debug flag is set, retains only the coefficient of e°, i.e. the finite part. It 
also applies a (mostly empirical) blend of simplification functions to get the compact expression for 
code generation. 

Two of these simplifications merit a closer look, for they affect also the numerical precision of 
the final result. The difference and sum of the sfermion mass squares not only appear frequently in 
the amplitude, they are also computed as a by-product of diagonalization anyway, with effectively 
no loss in precision even for degenerate sfermion masses. We thus substitute 

MSf2 [3,f , 3 ] = MSf2[2 ,t,g\ - MSf2 [1 ,f , 3 ] , 

MSf2[4,f,3l = MSf2[l,f,3] + MSf2[2,t,fl']. 

Considering that in the gaugeless limit with vanishing bottom mass the remaining sbottom mass 
square is equal to the stop-sector breaking parameter, to? = to? , the difference of stop and sbottom 
mass squares can be computed with better precision, too, which is particularly relevant for the H ± 
self-energy where this difference appears in denominators. Hence we also introduce 


MSq2Diff [?’, j] = MSf2[*,4,3] - MSf2[j,3,3] . 


• Command line: scripts/6-comb se 
se as in Step 1. 

• Main packages used: packages/FinalSimp .m, packages/ExpandDel .m 

• Input: m/sE/4-simp.0, m/sE/4-simp.1, m/rc/* 

• Output: m/sE/6-comb 

4.8 Step 7: Generating Code 

The last step writes out optimized Fortran code using FormCalc’s code-generation functions [42]. 
Self-energies and tadpoles are arranged in three groups, such that FeynHiggs can also perform a 
partial evaluation depending on its flags (e.g. if 2 x 2 mixing only is requested): 

1. A0A0, HmHp, 

2. hOhO, HHHH, hOHH, hO, HH, 

3. hOAO, HHAO, AO. 

Abbreviations are introduced for both loop integrals and common subexpressions, which shortens 
the entire expression significantly. Some variables need to be disambiguated, for example the 
sfermion masses in the gaugeless limit must not use the same identifiers as other sfermion masses 
in FeynHiggs. 

Lastly, the handwritten part of the code from the ‘static’ subdirectory is copied to the Fortran 
output directory ‘f’, which is then complete and replaces the entire src/TwoLoop/TLsp branch of 
the FeynHiggs source tree. 

• Command line: scripts/7-code 

• Main packages used: FormCalc, packages/FCSettings ,m 

• Input: m/*/6-comb, static/* 

• Output: f/* 

5 Conclusions 

We have described the implementation of the O(a^) Higgs-mass corrections for the complex MSSM 
in FeynHiggs. The original calculation of Ref. [29] was restructured entirely for this purpose and 
the new code is included in the FeynHiggs source code, available at http://feynhiggs.de. The 
calculation is split into seven tasks which were explained in detail. While each step is specific to 
the corrections computed here the code can be used as a template due to its modular structure and 
adapted with little effort for similar calculations. 
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